function [ rz,uz,zl1n,zl2n,zlin,rr,ur,rl1n,rl2n,rlin ] = laxwendrof( T, n, v)
%用Lax-Wendrof格式解不同初值的线性气动方程并将。
%参数说明：
%   输入
%   T：计算时间。
%   n：单位长度上的网格数量。
%   v：网格比。
%   输出
%   rz,uz：周期初值计算结果，在区域{(x,t)|-1<x<1,0<t<1}内的结果。
%           rz,为rou的值，uz为u的值。
%   zl1n:周期函数的L1误差
%   zl2n:周期函数的L2误差
%   zlin:周期函数的Linf误差
%   rr,ur：Riemann初值计算结果，在区域{(x,t)|-1<x<1,0<t<1}内的结果。
%           rr,为rou的值，ur为u的值。
%   rl1n:Riemann数据的L1误差
%   rl2n:Riemann数据的L2误差
%   rlin:Riemann数据的Linf误差
h=1/n;
tau = v*h;
t = floor(T/tau);%所需时间步，ie迭代次数。
s = t;%对于非周期函数的时候其中一个方向需要的额外网格。

%先计算Riemann数据
length = 2*s+2*n;
rr = zeros(2*n,t+1);
ur = zeros(2*n,t+1);

rr(1:n,1)=1;rr(n+1:2*n,1)=3;
ur(1:n,1)=5;ur(n+1:2*n,1)=1;

w1 = zeros(length,t+1);
w2 = zeros(length,t+1);

w1(1:s+n,1)=3;w1(s+n+1:length,1)=2;
w2(1:s+n,1)=2;w2(s+n+1:length,1)=-1;

pre =1;now=2;
for p = 2:t+1
    pre = p-1;
    now = p;    
    for i = p-1:length-p+1
        a1(i) = (w1(i+1,pre)+w1(i,pre))/2;
        a2(i) = (w2(i+1,pre)+w2(i,pre))/2;        
    end
    for i = p:length-p+1
        temp = (a1(i)*(w1(i+1,pre)-w1(i,pre))-a1(i-1)*(w1(i,pre)-w1(i-1,pre)))*(tau^2/(2*h^2));
        w1(i,now)=w1(i,pre)-(w1(i+1,pre)-w1(i-1,pre))*(tau/(2*h))+temp;        
        temp = (a2(i)*(-w2(i+1,pre)+w2(i,pre))-a2(i-1)*(-w2(i,pre)+w2(i-1,pre)))*(tau^2/(2*h^2));
        w2(i,now)=w2(i,pre)+(w2(i+1,pre)-w2(i-1,pre))*(tau/(2*h))+temp;        
    end    
end
rr=w1(s+1:length-s,:)-w2(s+1:length-s,:);
ur=w1(s+1:length-s,:)+w2(s+1:length-s,:);

%Riemann数据的真解：
nn = 2*n;
rrR = zeros(nn,t+1);
urR = zeros(nn,t+1);
for i=1:nn
    for j=1:t+1
        x = (i-n+1/2)*h;
        tt = tau*j-tau;
        if x<=-tt
            rrR(i,j) = 1;
        else if x <= tt
                rrR(i,j) = 4;
            else
                rrR(i,j) = 3;
            end
        end            
        if x<=-tt
            urR(i,j) = 5;
        else if x <= tt
                urR(i,j) = 2;
            else
                urR(i,j) = 1;
            end
        end            
    end
end
err = rr-rrR;
%mesh(err)
%pause
Size =size(err);
S = Size(1)*Size(2);
rl1n = sum(sum(abs(err)))/S
rl2n = sqrt(sum(sum((err).^2))/S)
rlin=max(max(err))

%计算周期函数
rz = zeros(2*n,t+1);
uz = zeros(2*n,t+1);
w1 = zeros(2*n,t+1);
w2 = zeros(2*n,t+1);

for i =1:nn
    x=(i-n-1/2)*h;
    w1(i,1)=(4+cos(pi*x)+sin(pi*x))/2;
    w2(i,1)=(cos(pi*x)-sin(pi*x))/2;     
end

pre =1;now=2;
for p = 2:t+1
    pre = p-1;
    now = p;
    for i = 2:nn-1
        w1(i,now)=w1(i,pre)-(tau/h)*(w1(i,pre)-w1(i-1,pre));
        w2(i,now)=w2(i,pre)+(tau/h)*(w2(i+1,pre)-w2(i,pre));
    end    
        w1(1,now)=w1(1,pre)-(tau/h)*(w1(1,pre)-w1(nn,pre));    
        w1(nn,now)=w1(nn,pre)-(tau/h)*(w1(nn,pre)-w1(nn-1,pre));    
        w2(1,now)=w2(1,pre)+(tau/h)*(w2(2,pre)-w2(1,pre));
        w2(nn,now)=w2(nn,pre)+(tau/h)*(w2(1,pre)-w2(nn,pre));
        %plot(w1(:,now));
        %pause(0.1);
end
rz=w1-w2;
uz=w1+w2;

%周期函数的真解：
nn = 2*n;
rzR = zeros(nn,t+1);
uzR = zeros(nn,t+1);
for i=1:nn
    for j=1:t+1
        x = (i-n-1/2)*h;
        tt = tau*j-tau;
        rzR(i,j) = 2 + cos(pi*tt)*sin(pi*x)+sin(pi*tt)*sin(pi*x);
        uzR(i,j) = 2 + cos(pi*tt)*cos(pi*x)-sin(pi*tt)*cos(pi*x);
    end
end
err = rz-rzR;
%mesh(err)
%pause
Size =size(err);
S = Size(1)*Size(2);
zl1n = sum(sum(abs(err)))/S;
zl2n = sqrt(sum(sum((err).^2))/S);
zlin=max(max(err));



end